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The effect of large scales on the statistics and dynamics of turbulent fluctuations is studied 
using data from high resolution direct numerical simulations. Three different kinds of forcing, 
and spatial resolutions ranging from 256 3 to 1024 3 , are being used. The study is carried out 
by investigating the nonlinear triadic interactions in Fourier space, transfer functions, structure 
functions, and probability density functions. Our results show that the large scale flow plays an 
important role in the development and the statistical properties of the small scale turbulence. The 
role of helicity is also investigated. We discuss the link between these findings and intermittency, 
deviations from universality, and possible origins of the bottleneck effect. Finally, we briefly describe 
the consequences of our results for the subgrid modeling of turbulent flows. 



I. INTRODUCTION 

The nature of the interactions in a turbulent flow has 
been a long standing problem. Most investigations are fo- 
cused on statistically homogeneous and isotropic turbu- 
lence. However, in most cases in nature this assumption 
is not necessarily true. Instead, one usually finds turbu- 
lence to be embedded in a large scale flow. In these cases, 
turbulence originates from some instability (e.g. convec- 
tion) and is not uniformly distributed in space and time. 
To what extend the flow is approaching a statistically 
homogeneous and isotropic state as smaller and smaller 
scales are developed is still an open problem. 

Due to insufficient computational resources (e.g in as- 
trophysics, weather and climate prediction, for tokamaks, 
and in industrial applications) the computational efforts 
to study turbulent flows at these high Reynolds numbers 
are compelled to resort to modeling of the smaller scale 
turbulent fluctuations for which numerical resolution is 
not available. However, in order to accurately model the 
small scales a good understanding of the impact of the 
large scale flow in the small scales is needed, in addition 
to the knowledge of what (if any) is the feedback of the 
small scale turbulent flow on the large scale structures. 
This latter process is usually modeled by an effective vis- 
cosity, or other forms of subgrid modeling of the Reynolds 
stress tensor. These two types of interactions (from large 
to small and from small to large scales, or "downscal- 
ing" and "upscaling" ) are often considered to take place 
only through a local cascade of energy, although there is 
evidence that nonlocal interactions between widely sep- 
arated scales are also relevant. 

Studies of local and nonlocal interactions have been 
done previously in direct numerical simulations (DNS), 
although at moderate spatial resolutions and Reynolds 

numbers 0,EHSIlliHHH[li[l2. WhUe S ° mG 
of these studies supported the existence of nonlocal in- 
teractions, it was argued by later studies that their pres- 
ence was associated with moderate values of the Reynolds 
number in the simulations, or that it was linked to the 
precise definitions used for transfer functions and the in- 
terpretation of the results. In some cases [E IS large 
eddy simulations (LES) were also used to extend the 



range of Reynolds numbers studied, although if nonlocal 
interactions are present the impact of the subgrid model 
on the transfer is unclear and is a point of study in it- 
self. Some of these studies also considered the effect of 
anisotropic or coherent large scale forcing (see e.g. Refs. 
Dl EJ El)- Evidence of nonlocal interactions between 
large a nd small scales has been found also in experiments 
|lall3lll4j . Observations of the persistence of anisotropy 
at small scales in experiments at large Reynolds numbers 
|l5l [T^ | , in the atmosphere |l7j , and in numerical simu- 
lations (see e.g. Ref. 01 an d references therein) also 
suggest the existence of a direct coupling between dis- 
parate scales. The presence of nonlocal interactions has 
also been associated in numerical studies with departures 
from universality and the development of intermittency 
[l9j . Recently, a new study at high resolution [2(| pre- 
sented a detailed analysis of nonlocal interactions using 
DNS for large Reynolds numbers, although for a partic- 
ular and coherent forcing function. 



In this work we examine three different kinds of flows 
in DNS in triple periodic boxes, generated by different 
body forces. We explore external forcings with and with- 
out helicity, and forcings with either infinite or short cor- 
relation times. The energy injection scale is varied, as 
well as the viscosity and spatial resolution of the runs, to 
explore a wide range of Reynolds numbers and configu- 
rations. In all cases, we observe the presence of nonlocal 
interactions coexisting with a local direct cascade of en- 
ergy. The intensity of the nonlocal interactions depends 
on the transfer function studied. For a particular (highly 
symmetric) forcing, we also observe a correlation between 
regions of large scale shear and strong small scale gradi- 
ents, supporting previous studies that linked intermittent 
effects with interactions with the large scale flow. 
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II. SETUP AND THEORY 
A. Equations, code, and simulations 

For an incompressible fluid with constant mass density, 
the Navier-Stokes equations are: 

d t u + u • Vu = -Vp + ^V 2 u + f , (1) 



V • u = 0, (2) 

where u is the velocity field, p is the pressure divided by 
the mass density, and v is the kinematic viscosity. Here, f 
is an external force that drives the turbulence. The mode 
with the largest wavevector in the Fourier transform of 
f is going to be denoted as kp and we are going to refer 
to 2-Kkp 1 as the forcing scale. We also define the viscous 
dissipation wavenumber as k v = (e/V 3 ) 1 / 4 , where e is the 
energy injection rate (as a result, the Kolmogorov scale 
is 7] = 2-7T I k v ) . A large separation between the two scales 
(kp 1 ;§> fc^ 1 ) is required for the flow to reach a turbulent 
state. 

In the absence of external forcing and viscosity, the 
Navier-Stokes equations in three dimensions have two 
ideal invariants: the energy 

E=\j u 2 d^ , (3) 

and the helicity 

H= i J u- V x udx 3 . (4) 

All the results discussed in the following sections result 
from analysis of data from direct numerical simulations 
of the Navier-Stokes equations. We solve Eqs. and 
(J2J) using a parallel pseudospectral code in a three dimen- 
sional domain of size 2tt with periodic boundary condi- 
tions [21], The pressure is obtained by taking the 
divergence of Eq. l(T]l. using the incomprcssibility con- 
dition (J2J), and solving the resulting Poisson equation. 
The equations are evolved in time using a second order 
Runge-Kutta method. The code uses the 2/3-rule for 
dealiasing, and as a result the maximum wavenumber is 
kmax = N/3 where N is the number of grid points in each 
direction. All simulations presented arc well resolved, in 
the sense that the dissipation wavenumber k v is smaller 
than the maximum wavenumber k max at all times. 

The Reynolds number is defined as R e = UL/v, where 
U is the r.m.s. velocity and L is the integral lengthscale 
of the flow 

f E{k)k~ l dk 
L = 2 - fE(k )d k > ® 

where E(k) is the energy spectrum. The large scale 
turnover time can then be defined as T = U/L. We 



can also introduce the Taylor based Reynolds number 
R\ = UX/v, where the Taylor lengthscale A is given by 

Several simulations were done with different resolu- 
tions (from A" = 256 to 1024) and kinematic viscosities. 
Table [I] shows the parameters for all the runs. The rms 
velocity in the steady state of all the runs is U ~ 1. To 
asses the effect of different large scale stirring forces, we 
used three expressions for the external force f . The first 
expression corresponds to a Taylor-Green (TG) flow [23l | 

frG = fo [sin(kFx) cos(k F y) cos(k F z)x— 

— cos(k F x) sin(k F y) cos(k F z)y] , (7) 

where fo is the force amplitude. This expression is not 
a solution of the Euler's equations, and as a result small 
scales are generated fast when the fluid is stirred with 
this forcing. The resulting flow has no net helicity, al- 
though regions with strong positive and negative helicity 
develop. 

In order to study directly the effect of helicity and its 
transfer between different scales, wc also did simulations 
using the Arn'old-Childrcss-Bcltrami (ABC) forcing 

f*ABC = fo{[Bcos(k F y) + Csin(k F z)}x+ 
+ [A sin(kpx) + C cos(kpz)] y + 
+ [A cos(k F x) + B sin(fc F y)] z} , (8) 

with A = 0.9, B = 1, and C = 1.1. The ABC flow 
is an eigenfunction of the curl and an exact solution of 
the Euler equations. When the flow is forced using this 
expression, turbulence develops after an instability sets in 
|24j . Unlike TG forcing, ABC forcing injects net helicity 
into the flow. 

The amplitude and phase of these two forcings is kept 
constant during the simulations, and as a result the ex- 
ternal force has an infinite correlation time. It is a com- 
mon practice in studies of isotropic and homogeneous 
turbulence to force in Fourier space injecting energy in 
all modes in a Fourier shell and changing the phase of 
each mode with a short correlation time. To compare 
with the results of TG and ABC forcing, we also imple- 
mented a random forcing 

fRND = / * X « i(k ' X+W (9) 

k|— fc_p 

where the phase 4>k of each mode with wavevector k was 
changed randomly with a correlation time r c that was 
taken to be t c = 0.1T. 

B. Scale interactions and transfer 

To investigate the interactions between different scales 
we split the velocity field into spherical shells in Fourier 
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TABLE I: Parameters used in the simulations. N is the linear 
grid resolution, f the forcing [either Taylor Green (TG), Bel- 
trami (ABC) or random (RND)], ItF the forcing wavenumber, 
v the kinematic viscosity, R e the Reynolds number, and R\ 
the Taylor based Reynolds number. 



Run N 



I 



256 TG 



II 512 

III 1024 

IV 256 

V 256 

VI 512 
VII 



TG 
TG 
ABC 
ABC 
ABC 



1024 ABC 
VIII 256 RND 



~2 2 x ltr 3 

2 1.5 x 10 -3 

2 3 x 10~ 4 
10 2.5 x 10~ 3 

3 2 x 10" 3 
3 6.2 x 10~ 4 
3 2.5 x 10" 4 
1 1.5 x 10~ 3 



Re R\ 



675 300 

875 350 

3950 800 

275 230 

820 360 

2520 670 

6200 1100 

2030 650 



space of unit width, i.e. v = Y^k w k wnere v k is a 
filtered velocity field such that only the Fourier modes in 
the shell K < \k\ < K + 1 (from now on called shell K) 
are kept. From equation the rate of energy transfer 
T%{K, P, Q) (a third-order correlator) from energy in shell 
Q to energy in shell K due to the interaction with the 
velocity field in shell P is defined as usual pol I25I [26ll2^ | 



T 3 (K, P,Q) = -J v K • (v P • V)v Q dx 3 



(10) 



Note that this term does not give information about the 
energy the shell P receives or gives to the shells K and 
Q. The computation of T3 for the three shells K, P, and 
Q from 1 up to a wavenumber K max requires ~ K^^N 3 
operations and is therefore demanding in computer re- 
sources. For example, in the 1024 3 simulations, to com- 
pute T^{K, P, Q) for all wavenumbers up to K max = 80, 
it takes as much computing time as evolving the hydro- 
dynamic code for two turnover times. 

If we sum over the middle wave number P we obtain 
the total energy transfer T^-ftT, Q) from shell Q to shell 
K: 

T 2 (K, Q)=J2 P, Q) = - j w K ■ (v ■ V)v Q dx 3 . 

p ^ 

(11) 

Positive transfer implies that energy is transfered from 
shell Q to K, and negative from K to Q; thus, both T3 
and To are antisymmetric in their (K, Q) arguments (see 
Ref. [27j). T 2 {K,Q) gives information on the shell-to- 
shell energy transfer between K and Q, but not about 
the amplitude of the triadic interactions themselves. 

If we further sum over the wave number Q we obtain 
the transfer 

Ti (K) = T 2 (K, Q) = -f vk ■ (v • V)vdx 3 (12) 



that expresses the rate the shell K receives energy from 
the velocity field (all shells). 



The energy flux is reobtained from these transfer func- 
tions as 

n ( fc ) = - Tl ^) = - / v a-( v • v)vrf (13) 



K=0 

where the notation 

K 



= Y v A " and = ^ v K , (14) 



K'=0 



K'=K+1 



is used (see Ref. |2Sj). 

We can further define the flux of energy at some given 
scale due to the interactions with some other scale as 

Hp(fc) = -J2 E T 3(^' P >Q) = - / v£-(v P .V)v>dx 3 

(15) 

which expresses the flux of energy at the scale K 1 due 
to the interactions with the scale P" 1 or in other words 
the energy flux at the scale K ~ 4 if the velocity field was 
advected just by the velocity field at scale P -1 . As will 
be discussed later, the question of locality depends on 
which of the different transfer functions or flux is under 
investigation. 

It is worth noticing that through the study of the three 
functions Ti, T2, T3, and the fluxes lip, we can measure 
the nature and intensity of the interactions directly, and 
therefore we avoid introducing a non-locality parameter 
as was done e.g. in Refs. 

The total energy balance equation for a given shell is 
written as 

d t E{K) = Ti{K) + vV{K) + F(K) (16) 
where we have also introduced the dissipation function 



vV{K) = vj |Vu A -| 2 dx 3 



(17) 



and the energy injection rate to the velocity field through 
the forcing term 



T(K) 



f • u K cfac 3 



(18) 



Finally, we will also investigate the transfer of helicity 
among different scales. Taking the inner product of the 
Navier-Stokes equation 0) with the vorticity at a scale 
K^ 1 (wk = V x vljc) and adding the inner product of 
the velocity at the same scale K -1 with the curl of (JTJ 
and space averaging we obtain: 

d t H{K) = Y.Q I ™k ■ (u x w Q )dx 3 + 

/ wk • Fdx. 3 +uj w K ■ V x w K dx. 3 . (19) 

Each term of the sum in the first term in the r.h.s of the 
equation above can be written as: 



T H (K, Q) = J • (u x w Q )dx 3 



(20) 




FIG. 1: (Color online) Rendering of enstrophy density in a 
1024 x 64 x 1024 slice of Run III. Bands of strong (1,2,3,4) 
and weak (A,B,C,D) shear in the external Taylor-Green force 
are indicated. The bars at the bottom show respectively the 
integral, Taylor, and dissipation scales. 

Note that the anti-symmetric property Th(K,Q) = 
-T H (Q,K) holds for T H (K,Q) i.e. the rate the shell K 
is gaining helicity from the interaction with the field wq 
and u is equal to the rate the shell Q is loosing helicity 
through the same interaction. This allows us to interpret 
Th(K,Q) as the transfer of helicity from the scale Q _1 
to the scale K~ x (see j2t| for the equivalent definition of 
the magnetic helicity transfer). 

III. INFLUENCE OF THE LARGE SCALE 
FLOW ON TURBULENCE 

We discuss in this section properties of Run III. First 
we describe the geometry of the resulting flow and we 
present general statistical results, such as probability 
density functions and structure functions. Then, we ana- 
lyze in detail the scale interactions using the formulation 
discussed in the previous section. Preliminary results for 
this flow were presented in [2(j. The results discussed 
here will be used as a reference to compare with the rest 
of the simulations in the following sections. 



A. Statistical properties, structure functions, and 
scaling exponents 

Run III is a 1024 3 simulation using TG forcing. Af- 
ter reaching a turbulent steady state, the simulation was 
continued for 10 turnover times. The expression of the 



4 




FIG. 2: (Color online) Zoom on a region with large enstrophy 
density showing field lines inside (above) and in its surround- 
ing (below) in Run III. Velocity field lines are helical inside 
and in the vicinity of the tube. Note the merging of two 
tubes in the south west corner. The side of the region shown 
is approximately one tenth of the side of the total domain. 



forcing has several spatial symmetries, and some of these 
symmetries are recovered in the flow in a statistical sense. 
Figurc^shows a rendering of enstrophy density in a thin 
slice of 1024 x 64 x 1024 grid points. Thin and elongated 
structures (vortex tubes) can be identified. As observed 
in previous studies, these structures are characterized by 
one large lengthscale in one direction (ranging from the 
Taylor to the integral scale) and a short lengthscale (close 
to the viscous dissipation scale) in the other two direc- 
tions. 

Since kp = 2, there are four planes in the z direction 
where the external force is zero (see Eq. [7J and the shear 
is maximum. In the steady state turbulent flow, these 
planes can be easily recognized since four bands where 
the vorticity is stronger are formed around them. Regions 
with less (or weaker) vortex tubes separate these planes, 
regions centered around the planes where the large scale 
shear has a minimum. Together, these two sets of regions 
form a large scale pattern that is observed to persist for 
long times, giving four "quiet" stripes and four stripes 
of "stronger" turbulence (the boundary between these 
regions is not well defined and fluctuates in space and 
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FIG. 3: Energy spectrum compensated by k 5 / 3 for Run III. 
The inset shows the energy flux. 

time, see Fig. QJ. Considering the persistence of these 
statistical symmetries of the flow, the velocity increments 
for this Run will be computed for displacements only in 
the x-y plane. 

When individual vortex tubes are studied, it is seen 
that the flow inside and surrounding the vortex tube is 
helical, as noted before by several authors [3(], El HE 
IsH l34j . Figure shows a rendering of a small region 
in the domain (approximately one tenth of the box). A 
region of large enstrophy density is indicated by dark 
(blue) colors, and field lines inside and in the vicinity 
of the enstrophy containing region are shown. In this 
region and its surroundings, field lines are helical, while 
the flow far from the region is not. This feature has been 
verified for a large set of vortex tubes in the complete 
domain. It has been claimed in the past [ISOmI that the 
development of these helical structures in a turbulent flow 
can lead to the depletion of nonlinearity and a quenching 
of local interactions. We will come back to this issue in 
the following sections. 

In spite of the presence of the large scale pattern, many 
features often associated with isotropic and homogeneous 
turbulence can be observed in this simulation. Figure |3| 
shows the angle-averaged Fourier energy spectrum and 
energy flux. An inertial range with constant energy flux 
is observed, together with a Kolmogorov-like scaling and 
a bottleneck as the dissipative range is reached. When 
probability density functions (pdfs) of transverse velocity 
increments 

6v±(x.,l) = r x [v(x + Zr) - v(x)] , (21) 

are computed in the whole domain (Fig. [3J, we observe 
distributions close to Gaussian for large increments, and 
the development of exponential and stretched exponen- 
tial tails as the increment I is decreased. 

Figures |S] and El show pdfs of transverse velocity incre- 
ments with I = 10?7 and 2rj respectively, but discriminat- 
ing between bands of strong and weak shear (as indicated 
in Fig. Bands with strong shear (1,2,3, and 4) have 
slightly but systematically stronger tails (i.e. a larger 




FIG. 4: Probability density functions of transverse velocity 
increments for Run III, for increments I — 2O77 (•••)) ' = IO77 

( ), I = 477 ( ),l = 2t}( ), and I = n (solid 

line), where n is the Kolmogorov dissipation scale. 
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FIG. 5: Probability density functions of transverse velocity in- 
crements with I — Wn for Run III. The solid line corresponds 
to the whole domain, dashed lines to regions 1 to 4 (strong 
shear), and dotted lines to regions A to D (weak shear; see 
Fig. 0. Notice the weak but systematic differences: the four 
dashed lines (some of them overlapping) have slower decaying 
tails and the dotted lines have faster decaying tails. 



probability of strong gradients) than regions A,B,C, and 
D. Note this behavior is true for each individual region, 
and the difference increases as the increment I is de- 
creased. This confirms statistically what can be inferred 
from Fig. ^ regions of strong shear have a larger density 
of vortex tubes and stronger gradients. 
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<5v(l = 2?7) 

FIG. 6: Probability density functions of transverse velocity 
increments with / = 2-q for Run III. Same labels as in Fig. 
|S| Note that the differences between quiet and strong regions 
appear even more systematic at this scale than in Fig. 
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FIG. 7: First order structure function Si (I) for Run III. The 
solid line corresponds to the whole domain, dashed lines to 
regions 1,2,3, and 4, and dotted lines to regions A,B,C, and D. 
The Kolmogorov scaling and the linear slope corresponding 
to a smooth flow are shown as a reference. The arrow indi- 
cates the Kolmogorov scale rj. Notice again the systematic 
differences between the two regimes identified in Fig. 1. 



We can also compute longitudinal velocity increments 



Sv\\ (x, = f • [v(x + If ) - v(x)] . (22) 

Figures and |H1 show respectively the second and third 
order longitudinal structure functions, where the struc- 
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FIG. 8: Second order structure function £>2(0 for Run III. 
Labels are as in Fig. |7| 

ture function of order p is defined as 

£ p (0 = (*«J|(x,0). (23) 

At scales smaller than the dissipation scale the field is 
smooth and the structure function of order p scales as l p . 
Kolmogorov's 1941 theory of turbulence predicts in the 
inertial range a scaling S p (l) ~ P^ 3 , although corrections 
due to intermittency are known and in general S p (l) ~ 
fa, where ( p ^ p/3 are the scaling exponents. 

As with pdfs of velocity increments, a clear trend sep- 
arating regions of strong and weak shear is observed. In 
the range of scales corresponding to the inertial range, 
the four regions with strong shear (regions 1,2,3 and 
4) show a larger slope than the four regions with weak 
shear. The slope of the structure function computed in 
the whole box lies between these two values. Consider- 
ing the results from both pdfs and structure functions, 
a correlation is observed between small scale gradients 
and large scale shear. Note that a correlation between 
stronger tails in the pdfs of velocity increments and vor- 
tex tubes has already been observed in E3- Here, 
a correlation between these quantities with large scale 
shear is further observed, in agreement with Ref. [l^ 
that suggests intermittency is related with interactions 
with the large scale flow. 

B. Interactions, energy transfer, and flux 

FigureElshows the functions T 3 (K, P, Q) and T 2 (K, Q) 
evaluated at Q = 40 in the turbulent steady state of Run 
III. When triadic interactions between shells are stud- 
ied, the strongest interactions are with the shell P = 3, 
where the large scale forcing is. Local interactions be- 
tween Fourier shells (K ~ P ~ Q) are two orders of mag- 
nitude smaller than nonlocal interactions with P ~ hp- 
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FIG. 9: (Color online) Above: amplitude of triadic interac- 
tions as given by T 3 (K,P,Q) at Q = 40 for Run III. Solid 
lines correspond to positive contour levels, and dotted lines 
to negative values. The dark (red) shaded region indicates 
shells in P close to the forcing shell, while the light (orange) 
shaded region indicates shells in P outside the local octave. 
Below: total shell-to-shell energy transfer T2(K,Q), trans- 
fer T2 S (K, Q) due to interactions with P close to the forcing 
shell (dashed line) , and difference between these two transfers 
(dash-dotted line). 

As a result, when T^(K, P,Q) is summed over all P to 
obtain the shell-to-shell energy transfer T 2 (K, Q), the en- 
ergy is observed to be transfered to small scales locally 
between shells, but with a step proportional to kp. In 
T 2 (K, Q), the negative peak at K ~ Q — Uf means energy 
is transfered from this shell to the shell Q, while the pos- 
itive peak at K ~ Q + kp indicates energy is transfered 
to this shell from the shell Q. 

The function Ts(K, P,Q) was studied in the steady 
state of this run also for Q = 10 and Q = 20, obtaining 
the same quantitative results: a dominance (when com- 
pared in amplitude) of triadic interactions with the large 
scale flow over local triadic interactions. The T2(K,Q) 
function for all values of Q between 10 and 80 also peaks 
at K fn Q ± kp (sec Ref. |20|). These peaks at fixed 
values of K — Q in T 2 (K, Q) are just the signature of the 
strong triadic interactions with the large scale forcing. If 
we compute the shell-to-shell transfer between shells K 
and Q due only to interactions with the large scale flow 

6 

T 2 LS {K,Q) = Y,T 3 (K, PQ), (24) 
p=o 

we obtain most of the shcll-to-shcll transfer (see Fig. 
The remaining transfer T 2 {K, Q) - T 2 LS (K, Q) still peaks 
at larger wavenumbers, although it still does not peak at 
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FIG. 10: Total energy flux U(k) in Run III (solid line), flux 
n LS (fc) due to interactions with the large scale flow (dashed 
line), and nonlocal flux U NL (k). 



K ~ Q = 40. This is due to nonlocal interactions outside 
the octave band not related with the large scale forcing 
(light shaded region in Fig. EJ. Note that the definition 
of the range summed over P in Eq. I|24|) . or of octave 
bands (here defined as Q/2 < K < 2Q) is somewhat 
arbitrary. 

As previously mentioned, while individual triadic in- 
teractions as described by T3 (K, P, Q) arc dominantly 
nonlocal, the shell-to-shell transfer T 2 (K, P, Q) describes 
a local transfer of energy although through interactions 
with the large scale flow. Clearly when we sum over 
shells, the larger number of modes in the small scales 
start dominating over the large scale modes. Summing 
further over K and Q, from Eq. (|13|) we can obtain the 
total energy flux. In analogy with Eq. (|24|) we can also 
define the flux due to interactions with the large scale 
flow (dark shaded region in Fig. 

6 

n LS (fc) = ]Tnp(fc), (25) 

P=0 

and the nonlocal flux due to interactions outside an oc- 
tave band (light shaded regions in Fig. |5J) 

k/2 

n NL (fc) = £n P (A0, (26) 

where ITp(fc) is defined in Eq. (15). 

Figure [Till shows these three fluxes in the steady state 
of Run III. The large scale flow is only responsible for a 
small (but not insignificant) fraction of the total energy 
flux (H LS ~ 0.2II). Furthermore although the ampli- 
tude of the nonlocal flux II NL (fc) depends on the defi- 
nition of an octave band, it is remarkable that it peaks 
at wavenumbers close to the peak of the bottleneck in 
the energy spectrum (see Fig. |3J- This gives a di- 
rect confirmation that the bottleneck is due to the de- 
pletion of the energy transfer due to local interactions 
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with K ~ P ~ Q. These interactions are inhibited be- 
cause of the presence of a numerical and/or viscous cut- 
off in wavenumber, and nonlocal triads become dominant 
SI HIE! EI- Indeed, in Ref. H it was shown using 
the eddy-damped quasinormal Markovian (EDQNM) ap- 
proximation that the bottleneck disappears if nonlocal 
interactions are excluded from the computation. Note 
that the spurious drop of the fluxes n LS and II NL at 
k = 80 in Fig. QJJis because the functions T 3 (Q,P,Q) 
and T 2 (K, Q) were only computed up to this wavenum- 
ber. 

Therefore there seems to be a hierarchy concerning 
the importance of nonlocality when one investigates the 
transfer functions. At the most basic level when one in- 
vestigates the transfer function T3 the interactions with 
the large scale flow are the strongest ones. When aver- 
aged over the middle wave number the resulting energy 
transfer T2 becomes local but not self-similar with a max- 
imum at K — Q ~ kp. When averaged further to obtain 
the energy flux, the non local interactions with the large 
scale flow are only responsible for 20% of the flux. 



IV. SCALING WITH REYNOLDS 

In this section we discuss results from runs I, II, and 
III. The three runs are forced using Eq. J7J), and the 
only parameters changed between the runs are the spa- 
tial resolution and the kinematic viscosity. This inves- 
tigation allows us to study how the results about the 
non-local interactions change as we increase the Reynolds 
number. Qualitatively the functions T 3 (K,P,Q) and 
T2(K,Q) are similar for all three runs: the strongest in- 
teractions [Ts(K, P, Q)] are with the large scale flow and 
the energy transfer [Ts{K, Q)] is local with the two peaks 
at K— Q ~ kp. For this reason we focus here on the study 
of the local and nonlocal flux of energy as the Reynolds 
number is changed. 

In Figure ITT1 we show the ratio of the flux H hs (k) to 
the total flux II(fc). As discussed in the previous section, 
the flux due to interactions with the large flow gets di- 
minished as we sum over more and more modes. As a 
result, in the inertial range of Run III interactions with 
the large scale flow are only responsible for ~ 20% of 
the energy flux. In the smaller Reynolds number runs 
(I and II), the amount of flux due to interactions with 
the large scale flow is larger. This implies that as we 
increase the Reynolds number the fraction of the energy 
flux in the inertial range due to interactions with the 
large scale flow decreases. However, note that in Run 
III a region in the inertial range where the ratio n LS /IT 
is approximately constant is observed. This region with 
constant ratio n LS /IT is not present in the simulations 
with lower Reynolds number. Therefore DNS with even 
higher Reynolds number than what is accomplished here 
are needed to determine how does this ratio scale with 
the Reynolds number. 
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FIG. 11: Scale variation of the flux ratio II LS /II for increas- 
ing Reynolds numbers; Runs I (dash-dot line), II (dash line), 
and III (solid line) . Note the presence of a wavenumber range 
where II LS /II becomes constant only for the highest resolu- 
tion run. 
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FIG. 12: Compensated energy spectrum k 5 ^ 3 E(k) (solid 
line) and helicity spectrum k 5 ^ s H(k)/kF (dashed line) in the 
steady state of Run VII, where fc_F is the forcing wavenumber. 



V. EFFECTS OF DIFFERENT FORCING 
FUNCTIONS 

A. Forcing expression and correlation time 

Nonlocal interactions discussed in the previous section 
were observed for runs with large scale non-helical forc- 
ing with infinite correlation time. It is of interest to know 
how much of these results translate to other kind of forc- 
ings. In this section we compare results for runs III, VII 
and VIII. Run VII is a 1024 3 simulation with constant 
helical forcing, and as runs I-III, it also displays a well 
defined large scale flow (for a description of the ABC 
flow, see e.g. Ref. In Run VIII the phases of 

the forcing function are changed randomly with a short 
correlation time. Although the resolution in this run is 
smaller (256 3 ), we will focus here on the effect of the forc- 
ing, since a systematic study of changing the Reynolds 
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FIG. 14: Shell-to-shell energy transfer function Ti{K,Q) for 
Run VII. The superimposed curves correspond to different 
values of Q, with Q G [10, 80]. 



FIG. 13: Scaling exponents of velocity structure functions 
(p as a function of order p in Runs III (TG forcing) and VII 
(ABC forcing). The Kolmogorov K41 scaling 14911 is indicated, 
as well as the She-Leveque (SL) prediction l5fJl . Note the 
difference between TG and ABC. 



number for fixed forcing was presented in the previous 
section. 

Figure IT2"1 shows the energy and hclicity spectra com- 
pensated by a Kolmogorov law in the turbulent steady 
state of Run VII. As proposed in Ref. 1431 and observed 
in several simulations 0, 0, 0, l48j the spectrum 
of helicity follows a Kolmogorov law and is proportional 
to kpE(k). The transfer of hclicity will be discussed in 
Section ED 

Figure^Jshows the scaling exponents ( p of the longitu- 
dinal structure functions, computed for runs III and VII. 
The extended self-similarity (ESS) hypothesis 
was used to compute the anomalous exponents, which 
show similar behavior for both runs. However, the re- 
sulting exponents for p > 3 from Run VII (ABC forcing) 
are slightly smaller than the exponents from Run III (TG 
forcing); for example, in Run III £4 = 1.2737±0.0005 and 
C 8 = 2.136±0.006, while in Run VII C 4 = 1.2647±0.0005 
and Cs = 2.026 ±0.008. Note that the difference between 
the scaling exponents for TG and ABC at order p = 4, 
though small, is nevertheless more than order of magni- 
tude larger than the error arising from measuring them 
by virtue of the ESS hypothesis which leads to negligible 
errors. 

It should be noted that it is unclear whether this ob- 
served difference is due to problems linked with using 
the ESS methodology itself, or if it indicates a depar- 
ture from universality between these two flows because 
of loss of homogeneity, as exemplified in Fig. 1. Since 
the two forcing functions studied here are both to some 
degree anisotropic and inhomogeneous, subleading con- 
tributions due to departures from full symmetry could 
be responsible for this discrepancy. In this context, a de- 
composition of the structure functions into their isotropic 



and anisotropic components [Tsl IH3I IHil ] could help to 
asses the degree of universality of each component. These 
points will necessitate further study and better resolved 
flows. 

The shcll-to-shell energy transfer T 2 (K, Q) for Run VII 
is shown in Fig. |14l As for non-helical forcing (see 
Sec. IIIII and Ref. [27]]) the shell-to-shell transfer function 
peaks at K « Q ± kp, indicating the transfer of energy 
is local but mediated by nonlocal interactions with the 
large scales. If we sum over K and Q to obtain the en- 
ergy flux, we reobtain the results discussed in Sec. |^ 
for Run III: most of the energy flux is local, although 
in the inertial range ~ 20% of the total flux is due to 
interactions with the large scale forcing. 

The transfer function T2(K,Q) was also computed in 
the turbulent steady state of Run VIII (see Fig. H5la). 
In this run the phases of the external force are changed 
randomly with a short correlation time. It is noteworthy 
that even in this case with isotropic and random forcing, 
evidence is found of nonlocal interactions with the large 
forcing scale: the T2(K, Q) function peaks at K w Q±kp 
for all wavenumbers Q studied (note that in this run hp = 
1). Fig. El is discussed further below. 



B. Forcing scale 

Through all this work we have shown transfer functions 
indicating that the energy is locally transfered between 
scales but with a fixed step that we associated with the 
forcing scale. In this subsection we compare the trans- 
fer function T^JiT, Q) in simulations forced at different 
wavenumbers. 

Figure IT51 shows T2(K, Q) for several values of Q com- 
puted in the steady state of runs VIII (kp = 1), IV 
(hp = 2), and V (kp = 10). A clear correlation is ob- 
served between the position of the peaks (\K — Q\) in the 
shell-to-shell transfer function, and the forcing wavenum- 
ber kp. In the last case, the peaks at \K — Q\ slightly 
smaller than kp could be associated with low Reynolds 
number effects. 
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FIG. 17: Energy spectrum compensated by fc -5 / 3 for Runs V 
(dashed line), VI (dotted line), and VII (solid line). All runs 
are for ABC forcing at different Reynolds numbers. Note the 
wavenumber on the x-axis has been divided by the dissipation 
wavenumber to make all dissipation ranges coincide. A slope 
of 1/3 close to the bottleneck (corresponding to a power law 
of fc _4/3 in the energy spectrum) is shown as a reference. 



FIG. 15: Shell-to-shell energy transfer function Tz(K,Q) for 
runs with different forcing wavenumbers Icf as indicated by 
the vertical lines; Run VIII with Icf — 1 (a), Run V with 
kF = 3 (b), and Run IV with /cf = 10 (c). In each panel, 
the several curves correspond to different values of Q in the 
inertial and dissipative ranges. Note that for each run, the 
peak of energy transfer is centered close to kF- 
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FIG. 16: The spectrum of relative helicity H(k)/kE(k) com- 
pensated by fc -1 for run VII. Note that the region where the 
spectrum is flat indicates a power law of fc -1 in the relative 
helicity. At large wavenumber, a fc 1//2 slope is shown only as 
a reference. The increase at large k is indicative of an excess 
of relative helicity in the small scales. 



VI. HELICITY TRANSFER 

The transfer of helicity is readily studied in Run VII, 
since the external forcing injects maximum helicity and 
all scales are dominated by the same sign of helicity. We 
will make no attempt here to separate the different signs 
of helicity in the simulations (see e.g. Refs. |45l l46l l55j ) . 



Before discussing the transfer of helicity, it is of in- 
terest to study spectral properties of helical flows. As 
previously discussed, the spectrum of helicity follows an 
approximate fc -5 / 3 law (see Fig. 1120 . As a result, the 
relative helicity H ( k )/ E { k ) k follows in the inertial range 
a fc -1 slope |2a, |43, Eg (he. small scales are less heli- 
cal than large scales). It was predicted in Refs. |5rl l57j 
that the dissipation scale of helicity should be larger than 
the energy dissipation scale, giving as a result a drop in 
the spectrum of relative helicity faster than fc -1 for small 
scales. This argument would be in agreement with the 
idea that small scales slowly recover the mirror symme- 
try broken by the injection of helicity in the large scales. 
However, previous simulations |47j and this high resolu- 
tion run both suggest that there is an excess of relative 
helicity in the small scales, when compared with the A; -1 
drop. Note that this slower than predicted recovery of 
symmetries in the small scales is also in agreement with 
the slower than ex pec ted recovery of isotropy observed 
in experiments [TR llfij. 

Figure IT7)1 shows H(k)/E(k), i.e. the spectrum of rel- 
ative helicity R(k) = H{k)/kE{k) compensated by fc -1 . 
A scaling of fc -1 for the relative helicity R(k) thus corre- 
sponds to a flat spectrum in Fig. 16, as observed through 
the inertial range up to k ~ 20. However, at small scales 
the compensated spectrum of relative helicity grows, pos- 
sibly as fc 1 / 2 or steeper, indicating that the spectrum of 
helicity H(k) at small scales is dropping slower than the 
spectrum of energy E{k) (see also Fig. EJ. This can 
be associated with the presence of vortex tubes, that are 
known to be helical (see Fig. [2J|. 

The effect of helicity, and indirectly the effect of non- 
local interactions which can be modeled for example 
through the introduction of a second time-scale in the 
problem in order to distinguish between the eddy turn- 
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FIG. 18: Shell-to-shell helicity transfer T H (K, Q) in Run VII, 
for all values of Q from 10 to 80. 



over time and a helical characteristic time, has also been 
used to explain the development of the bottleneck ef- 
fect that occurs at the onset of the dissipative range, see 
e.g. [f^l ■ Ref. [5^ predicts a k~ 4 / 3 energy spectrum 
for the bottleneck, that we found compatible with our 
spectra for runs V-VII (Fig. I17f) . We also observe a 
A: -4 / 3 range close to the bottleneck in the Taylor-Green 
flow that has no net helicity. Note that the argument 
in Ref. |H3 is based on the presence of non zero he- 
licity locally. The transition from the inertial range to 
the bottleneck seems to be dependent on the Reynolds 
number. It is worth noting that while at resolutions of 
256 3 mostly a bottleneck is observed, in the 1024 3 runs a 
short Kolmogorov-like scaling is found before the bottle- 
neck takes place. The origin of the bottleneck based on 
the dominance of nonlocal interactions in the dissipative 
range has also been investigated jH, |3^, ^| as we 
discussed in Sec. IIIII and is independent of the presence 
of helicity. These arguments are not mutually exclusive, 
as local in space generation of helicity in the small scales 
can also quench local interactions between eddies of com- 
parable sizes. However, the prediction in Ref. |39j for the 
spectral shape of the bottleneck is in disagreement with 
the spectra obtained in all simulations here. These points 
deserve further study at higher resolution if one is to be 
able to distinguish between the different ranges that may 
be occurring. 

Finally, we discuss the shell-to-shell helicity transfer 
T H (K,Q) in Run VII (Fig. [TSJl. As for the transfer of 
energy, the helical transfer function peaks at K w Q±kp, 
although this transfer is noisier. The helicity cascades di- 
rectly to smaller scales as the energy, confirming previous 
studies 0IH, HI!!!]. Tne helicity is not a positive def- 
inite quantity, and as a result both signs of the transfer 
can indicate a direct cascade depending on the sign of 
the helicity itself. Moreover, as we study smaller scales 
and the relative helicity decreases, both signs have more 
predominance increasing the noise observed in Th(K, Q). 



VII. CONCLUSIONS 

We start by summarizing our results. The energy 
transfer, triadic interactions, and statistical flow prop- 
erties were studied in several high resolution DNS with 
periodic boundaries. The spatial resolution ranged from 
256 3 to 1024 3 . The Reynolds number based on the in- 
tegral scale spanned values from R e « 275 to 6200, 
while the Taylor-based Reynolds number varied between 
R\ « 230 and 1100. Three forcing functions were used, 
two coherent functions with zero and maximum net he- 
licity respectively, and a random forcing with a short cor- 
relation time. The energy injection scale was also varied 
to study its impact on the energy transfer. 

Most statistical studies were done for Run III, a high 
resolution simulation with R\ sa 800 and TG forcing. 
The forcing and the resulting flow have spatial symme- 
tries that allow us to identify planes with strong and weak 
large scale shear easily. While in the whole domain the 
standard results were rcobtaincd (e.g. exponential and 
stretched exponential tails in the pdfs of velocity incre- 
ments, and anomalous scaling of the structure functions), 
when studying individual regions a correlation between 
large scale shear and small scale gradients was found. In 
regions of strong shear, the tails of the pdfs of velocity in- 
crements are stronger than in regions of weak shear, even 
at scales as small as twice the dissipation scale. Also, 
structure functions show a slightly larger slope in regions 
of strong large scale shear, i.e. a larger departure from a 
Kolmogorov p/3 self-similar scaling. 

A correlation between stronger tails in the pdfs of ve- 
locity increments and the presence of vortex tubes has 
already been observed in different decompositions of the 
flow jH El- Here, we observed a correlation between 
these quantities and the large scale shear. The corre- 
lation was observed to be persistent even at scales as 
small as the dissipation scale. Small differences were also 
observed in the anomalous scaling of the structure func- 
tions for TG and ABC forcing. However, it is unclear 
whether this is related with non-universal effects asso- 
ciated with interactions with the large scale forcing, or 
with the use of the ESS hypothesis. For individual vor- 
tex tubes in several flows, we verified that the flow in- 
side and surrounding the vortex tube is helical, as found 
before HJ EJ EH HI HI . Note that the development 
of helical structures in a turbulent flow can lead to the 
depletion on nonlinearity and a quenching of local inter- 
actions [35ll36T|. 

Concerning the energy transfer and triadic interac- 
tions, we confirmed for several forcing functions that the 
cascade of energy is local between Fourier shells, although 
it is strongly mediated by individual triadic interactions 
which are nonlocal in nature. As a result, the energy 
cascades from one shell to the next with a fixed step pro- 
portional to the forcing wavenumber kp. This effect was 
observed even in simulations using random forcing with a 
correlation time one order of magnitude smaller than the 
large scale turnover time. No qualitative differences have 
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been observed as the Reynolds number was changed. 

However, the degree of nonlocality observed depends 
on the quantity studied. A hierarchy is found in the 
relative amplitude of nonlocal effects in triadic interac- 
tions, the shell-to-shell energy transfer, and the energy 
flux. Triadic interactions are dominated by interactions 
with the large scale how. As more modes are summed to 
define the shell-to-shell transfer and the flux, the larger 
population of small scale modes starts dominating, and 
in the 1024 3 simulations the large scale flow is responsible 
for only ~ 20% of the total flux. 

An increase of the relative contribution of nonlocal 
interactions (with the large scale flow and with modes 
outside the octave band studied) in the total flux is ob- 
served as the dissipative range is reached. This result is 
in good agreement with claims that the bottleneck effect 
is due to the quenching of local interactions at the end 
of the inertial range because of the presence of a cut-off 
in wavenumbers pa , I40L Elj . 

It is worth noting that the local energy cascade through 
nonlocal interactions, or non negligible nonlocal interac- 
tions, has been observed in the past in simulations at 
lower Reynolds numbers (see e.g. Rcfs. @,|HiJE3) Our 
results confirm the presence of interactions between dis- 
parate scales in a turbulent flow at much larger Reynolds 
numbers, and for a variety of forcing functions and forc- 
ing scales. The results discussed here also shed some light 
on the controversy in the literature about the relevance 
of the nonlocal interactions. It has been claimed that in- 
teractions are local if a different measure for the locality 
is introduced 0, , or if wavelets or a binning of Fourier 
space in octaves is used ^lj (note that wavelets naturally 
introduce a binning in octaves of the spectral space) . The 
hierarchy found in the different transfer functions is the 
reason for this apparent inconsistency between previous 
results. As more modes are summed (e.g. to define oc- 
taves in spectral space, or to define our partial fluxes 
Hp), the small scales overcome the triadic interactions 
with the large scale flow, and local interactions give the 
largest contribution. This is also in agreement with re- 
cent theoretical results about the locality of the energy 
flux [5^ , or the locality of the shell-to-shell energy trans- 
fer in Fourier octave bands [60L l6ll| . 

The fact that in simulations with R\ ~ 1000 most 
of the flux is due to local interactions, does not preclude 
however the existence of strong nonlocal interactions with 
the large scale flow at the triadic interaction level, and 
it should be kept in mind that even at large values of 
R\ these interactions are responsible for a non-negligible 
fraction of the total flux. The presence of nonlocal in- 
teractions are a deviation from the standard hypothesis 
often associated with Kolmogorov (K41) theory [43, and 
can possibly explain departures from self-similar models 
of the higher order structure functions and controversial 
results observed in experiments, such as the slower than 
predicted recovery of isotropy in the small scales . 

Similar results were obtained for the helicity transfer. 
The injection of net helicity in the large scales breaks 



down the mirror symmetry in the flow. While confirming 
the direct cascade of helicity [il lil HI lil lil l48j . we 

also found a slower than expected recovery of the symme- 
tries in the small scales, with an excess of relative helicity 
at small scales compared with the expected drop. As 
in the energy cascade, the cascade of helicity takes place 
in fixed steps proportional to the forcing scale, indicating 
strong nonlocal triadic interactions. 

Nonlocal interactions can also be responsible for ob- 
served departures from universality. In [2(J | it was shown 
with a simple dimensional argument how nonlocal inter- 
actions can still be consistent with a k~ 5 ^ 3 energy spec- 
trum. Vortex tube stretching by the large scale flow plays 
a significant role, and the argument has points in com- 
mon with multifractal models of intermittency, such as 
the (3 model (see e.g. 0)- As previously mentioned, in 
this work we showed evidence of a positive correlation of 
large scale shear and small scale gradients, both in pdfs 
of velocity increments and in structure functions. Ref. 
[l9j showed that the anomalous scaling of the structure 
functions is reduced when nonlocal interactions with the 
large scale flow are artificially suppressed in a simulation. 
Both results suggest that departures from K41 theory can 
be associated with the imprint of the large scale forcing 
in turbulence. 

The modern language used in turbulence is related to a 
great extent to the K41 theory developed for the isotropic 
and homogeneous case, with the assumption that such 
fundamental symmetries of the equations would be recov- 
ered in the small scales even when broken (e.g. through 
external forcing) in the large scales. To apply this the- 
ory to real turbulent flows, the idea that small scales 
restore isotropy and homogeneity is based strongly on 
the assumption of local interactions between scales. The 
persistence of anisotropics and other deviations from K41 
observed in experiments and simulations have been asso- 
ciated in this work with the presence of strong non-local 
triadic interactions. At this point, the reader could ask 
how much of the edifice of homogeneous and isotropic 
turbulence remains. 

One of the assumptions in K41 theory is that the prop- 
erties of the inertial range are universal. This is directly 
related to the assumption of local interactions. Since 
eddies in the inertial range only interact with eddies of 
similar size, as the energy cascades through the inertial 
range a self-similar solution is obtained. The fraction of 
the flux due to non-local interactions, found to be ~ 20% 
here, can be interpreted as a subleading contribution to 
the local flux. Deviations from K41 theory are well doc- 
umented and several theories have been proposed to ex- 
plain them (see e.g. |SE 121, EE E3, E2 ) They have also 
been derived theoretically for the passive scalar in the 
context of the so-called Kraichnan model |^. Several 
of these corrections were associated with viscous effects 
for turbulence at finite Reynolds number. What our re- 
sults show is that these deviations can come as well from 
interactions with the large scale flow. In this context, 
simulations at higher resolution could help to study the 
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scaling of the non-local flux with the Reynolds number. 

It is known that in the Kraichnan model only prefac- 
tors of the scaling laws depend on the large-scale forc- 
ing (i.e. the exponents are universal: independent of 
the forcing). However, in the case of the Navier-Stokes 
equations, this has not been proved yet since the hierar- 
chy of equations for the n-point correlation functions is 
not closed. Then, inertial range solutions with exponents 
that depend on the forcing may exist. 

In the case of realistic turbulent flows as encoun- 
tered in astrophysics and geophysics, the presence of 
strong non-local interactions can lead to the persistence 
of anisotropics in the small scales. In this case, the 
anisotropic effects question the applicability of the the- 
ory of isotropic and homogeneous turbulence to real flows 
(at least at Reynolds numbers comparable to the ones 
studied in this work). The properties of the large scale 
flow can be more important than expected to shape the 
small scales. A systematic study of the isotropic and 
anisotropic contributions to the scaling 0, 13 can 
be a first step to recognize universal (and non-universal) 
features in these cases. 

Finally, as previously indicated in [2(| , the existence of 
non-negligible nonlocal interactions in a variety of turbu- 
lent flows gives support to models involving as an essen- 
tial agent of the nonlinear energy transfer the distortion 



of turbulent eddies by a large-scale flow (e.g. as in rapid 
distortion theory and its variants or as in the La- 

grangian averaged Navier-Stokes equations where 
the turbulent flow interacts with a smooth velocity field). 
We believe that the understanding of scale interactions 
in turbulence can lead to the development of a new gen- 
eration of subgrid models, beyond the usual hypothesis 
of locality done in most Large Eddy Simulations. 
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